Boundary effects in the critical scaling of entanglement entropy in ID systems 
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We present exact diagonalization and density matrix renormalization group results for the entan- 
glement entropy of critical spin-1/2 XXZ chains. We find that open boundary conditions induce an 
alternating term in both the energy density and the entanglement entropy which are approximately 
proportional, decaying away from the boundary with a power-law. The power varies with anisotropy 
along the XXZ critical line and is corrected by a logarithmic factor, which we calculate analytically, 
at the isotropic point. A heuristic resonating valence bond explanation is suggested. 
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Introduction — Quantum spin chains have recently 
acquired the status of paradigmatic systems to investi- 
gate and understand the subtle interplay between quan- 
tum entanglement and quantum criticality 0, E3> El @ • 
While much theoretical work in this area has focussed on 
periodic boundary conditions (PBC), experimental sys- 
tems typically have open boundary conditions (OBC). 
In this letter, we study a new effect of OBC on quantum 
entanglement in antiferromagnetic chains. 

The entanglement entropy (EE), first introduced by 
Bennett et al. in a quantum information context, has 
recently received more attention from a condensed mat- 
ter point of view since it has been shown to exhibit a 
universal scaling for some one dimensional (ID) quan- 
tum critical points 0, 0] ■ Defined as the von Neu- 
mann entropy of the reduced density matrix p(x) of a 
given subsystem of size x embedded in a larger closed 
ring of size L, a universal scaling is expected at a quan- 
tum critical point in the limit x 3> 1 and L — > oo: 
S(x) = — Trpln/5 = | In a; + si. The prefactor c is the 
central charge of the associated conformal field theory 
(CFT) 0, and si is a non universal constant related 
to the ultra-violet cut-off 0. The divergence of S(x) 
with x reflects an interesting property of the reduced den- 
sity matrix: p(x) has a number of non negligible eigen- 
values which also diverges with the subsystem size, like 
M(x) ~ c 5 - x c / 3 . 

This critical scaling has been verified numerically 
by diagonalizing frce-fcrmion type Hamiltonians for the 
quantum Ising chain (c = 1/2) or the XX chain (c 
1 . An important extension to critical and non-critical 
systems with finite size, finite temperature and different 
boundary conditions has been achieved by Calabrese and 
Cardy |2|. Using CFT, they showed for instance that for 
critical systems of finite size L with OBC, the expression 
for the EE of a subsystem of size x including the open 
end becomes: 

S(x, L) = £ In fdn(™)) + hi .9 + fll /2, (1) 
where In g is boundary entropy introduced in 0- (We 



have corrected factors of 2 in the last 2 terms in Eq. 
following [ljj.) 

In the following, we investigate the EE for a large 
class of critical XXZ quantum spin chains of length L 
with OBC. In addition to confirming the behavior (I), 
at large x, we find an unexpected alternating term in 
the EE which decays away from the boundary with a 
power law. The boundary also introduces a slowly de- 
caying alternating term in the energy density. We show 
numerically that these are proportional. The microscopic 
origin of the even-odd alternation of the EE will be un- 
derstood through a qualitative resonating valence bond 
(RVB) picture. 

Critical chains with open ends — The Hamiltonian for 
XXZ quantum spin chains of length L with OBC is 



l-i r 



x=l 



-(S X S X+1 + S x S x+1 ) + AS X S X+1 



(2) 



The critical regime is achieved for |A| < 1. The energy 
density (h x ) = ((S+S x+1 + S x S+ +1 )/2 + AS^ +1 ) is 
independent of x in periodic chains. On the other hand, 
an open end breaks translational symmetry and the en- 
ergy density picks up a slowly decaying alternating term 
or "dimerization" : 



(h x ) = E u (x) + (-l)*E A (x), 



(3) 



where Ea(x) becomes non zero near the boundary and 
decays slowly away from it. We can calculate E\ (x) 
by Abelian bosonization methods modified by OBC jllj . 
The low energy effective Hamiltonian is simply that of 
a free masslcss relativistic boson, <j>{x) with a boundary 
condition, (f)(0) = constant. Eu(x) goes to a constant 
like 1/x 2 (for all |A| < 1) but this is not of interest here. 
The alternating part of h x is (— l) x+1 sin(V AttKc/)) where 
is the Luttinger liquid parameter. Us- 
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ing the standard conformal mapping we obtain in a finite 
system of length L 
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At the Heisenbcrg [SU(2) invariant] point, A = 1, the 
low energy effective Hamiltonian contains a marginally 
irrelevant coupling constant, A. It leads to a contribution 
to the scaling dimension of 0(A) for most operators 0, 

. The staggered energy density, sin y/2ir4> nas a scaling 
dimension lial 



7 = 1/2 + tt\/3A + 0(A 2 ). 



(5) 



This implies [l3| that the bulk correlation function G = 
(hihj) decays as l/[\i — j | (In \i — j\) 3 ^ 2 }- With a bound- 
ary, one can derive a renormalization group equation for 
E A {x): 

[d/dQnx) + /3(\)d/dX + j(X)]E A {x) = 0, (6) 

which is the same as that obeyed by G(x) except that 
the term 7 gets replaced by 27. (Since we are interested 
in Ea(x) for x ^> 1, the presence of the boundary does 
not affect the scaling dimension of sin(v27r^). This can 
be seen from considering the operator product expansion 
of the marginal operator with sm(y/2n(j>) . At short dis- 
tances, it is unaffected by the boundary. On the other 
hand, corrections to scaling dimensions of boundary op- 
erators due to the marginal operator are affected non- 
trivially by the boundary condition. See |l8j.) This im- 
plies Ea(x) oc l/[-\/[xf(ln |x|) 3//4 ]. It is highly non-trivial 
to include both the log corrections and the finite size ef- 
fects in Ea{x,L). However, there is a simple result at 
x = L/2. Including the cubic term in the /3-function for 
the marginal coupling constant, and other higher order 
corrections 0, this becomes: 



,L , 
e a{-^,L) = a 



l + a 2 /[\n{L/ ai )Y 



L[\n{L/ ai ) + (l/2)lnln(L/ai)p 



(7) 



where ao, ai and a 2 are constants. This formula is fit 
to density matrix renormalization group (DMRG) data, 
(after extracting uniform and alternating part by fitting 
both locally to a polynomial using a 5-point formula) in 
Fig. ^ obtaining excellent agreement. In the same figure 
we also show Ea(L/2) for the isotropic model with first 
and second neighbor interactions with J2/J1 = 0.241167. 
This model is at the critical point between gapless and 
gapped spontaneously dimerized phase and the marginal 
coupling constant, and hence the log corrections are ex- 
pected to vanish here 0], as is seen in the figure. 

RVB picture — This modulation in (h x ) can be inter- 
preted as an alternation of strong and weak bonds along 
the chain. Indeed, the boundary spin at x — 1 will have a 
strong tendency to form a singlet pair with its only part- 
ner on the right hand side. On the other hand the spin 
located at x = 2 will be consequently less entangled with 
its right partner at x = 3 since it already shares a strong 
entanglement with its left partner. This tendency to- 
wards dimcr formation over odd bonds induced by OBC 
can be naively depicted through the simple valence bond 
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DMRG m=512 
□ Jj/J, =0.241 167 

o J 2 /J n =o 



Eq. (7) (a =1.374, a, =0.377, a 2 =3.38) 




FIG. 1: E A (L/2,L) computed by DMRG (keeping m = 512 
states) at the SU(2) point up to L = 1000 for the near- 
est neighbor model J2/J1 — (o) and at the critical point 
J2/J1 = 0.241167 (□) (a). Full lines are fits, indicated on the 
main panel. Zooms on large L regions are showed in (b) and 
(c): (b) absence of log corrections for J2/J1 = 0.241167; (c) 
excellent agreement between Eq. (7) and DMRG. 
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FIG. 2: A particular valence bond state. The number of 
valence bonds crossing each link, N(x), is given. 
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t2x-l|2x} — |l2a;-lT2x)] 



(8) 



(This is the exact ground state of the Majumdar-Ghosh 
model with a second neighbor coupling obeying J 2 j J\ = 
1/2 0|.) In such a state, the EE of a finite subsystem 
alternates between In 2 if the interface between the sub- 
system and the rest of the chain cuts a singlet bond and 
otherwise. More generally, we may write any singlet state 
as a linear combination of all valence bond states (with 
no bonds crossing each other). It is trivial to calculate 
the EE for a simple state corresponding to any particu- 
lar valence bond configuration, such as the one drawn in 
Fig. (J2J. It is given exactly by S(x) = iV(x)ln2, where 
N(x) is the number of valence bonds crossing the link 
x between sites x and x + 1. Since one may think, in a 
renormalization group sense, of the ground state of the 
random-bond Heisenberg model as consisting of a sin- 
gle valence bond configuration 0, this leads 0] to an 
interesting prediction for the EE. Unfortunately, the non- 
random critical case, where the ground state involves an 
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infinite sum over different valence bond configurations, 
including bonds of arbitrary length, appears more diffi- 
cult. Nevertheless, it seems plausible that the enhanced 
tendency towards valence bonds on odd links induced by 
a boundary will translate into an alternating term in the 
EE. 

Numerical results for the entropy — It is well known 
that when the Ising exchange term A = 0, using a 
Jordan- Wigner mapping, the XXZ Hamiltonian (J2J) is 
equivalent to a free fermion Hamiltonian 



7~Lxx 



L-l 



x+l 



(9) 



which can be solved in momentum space. The computa- 
tion of the EE in this non-interacting case can be easily 
achieved numerically for very large systems since it only 
requires diagonalization of an L x L matrix (see Ref. |l6( 
for some details). We first report results from exact di- 
agonalizations of the free fermion Hamiltonian JHJ with 
L = 2000 sites. In Fig. we superimpose both PBC 
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FIG. 3: Results from exact diagonalization of XX chains with 
L = 2000 spins 1/2. EE 5(x,2000), plotted vs the subsystem 
size x in PBC (upper symbols) and OBC (lower symbols); 
the insets shows zooms around the chain center. The full 
lines are fits: Eq. JTOJ with si = 0.726067 for the PBC case 
whereas in the OBC case, some uniform and staggered terms 

0.055-(~l)^0.25 have been added tQ g «y 



and OBC cases. In the PBC situation, as predicted in 
Ref. , the entropy is very well described by the follow- 
ing expression 



S PBC (x,L) 
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(10) 



with c = 1 and Si ~ 0.726 [13] ■ On the other hand, 
the open chain does not obey the expression but 
as expected from the RVB picture presented above, an 



alternating term is found, as is visible in the lower inset 
of Fig. |3 Indeed, the strong odd bond - weak even bond 
picture agrees qualitatively with an enhanced (reduced) 
entropy for subsystems cutting an odd (even) bond. We 
thus define the uniform and alternating parts of the EE 
for large x: 



S(x,L)=Su(x,L) + (-l) x S A (x,L), 



(11) 



where both Sjj{x,L) and Sa{x,L) are expected to be 
slowly varying functions for x 3> 1. Interestingly, for the 
XX case, Sa(x, L) is found to decay slowly like a power- 
law, 



S^(x,L) 



1 



(12) 



which has the same exponent as the alternating part 
of the energy density Ea, as follows from Eq. (@J us- 
ing the correct value K = 1, for the XX model. It 
turns out that both the alternating part of the EE 
Sa and Ea, defined by Eq. Q, are nearly propor- 
tional at large x with a proportionality constant of 7r/2: 
S A {x,L) = (-ir/2)E A (x 1 L) + 0{l/x 2 ). S v (x,L) has also 
a boundary-induced correction to the result in Eq. 
which is also oc l/[Lsm(irx/L)]. This is included in the 
fit in Fig.El 



0.08 



0.06 



• 


A=-0.75 


♦ 


A=-0.5 


▼ 


A=0 


A 


A=0.5 


□ 


A=l 



< 

(/) 0.04 



0.02 



/4 



* P 



UJ 

C/3 




_F/4e-04 

/// s A 

///' 2e-04 



(a) 



-1 -0.5 0.5 

A 



_l_ 



0.02 0.04 0.06 0.08 



0.1 



0.12 



FIG. 4: Linear behavior of the alternating part of the EE, Sa 
as a function of the alternating energy density, — Ea, both 
computed using DMRG on critical open XXZ chains of size 
200 < L < 1000 for a few values of the anisotropy A. Data 
from ED at A = are also shown for L = 2000. Dashed 
lines are linear fits of the form Sa = —oEa (see text). The 
inset (a) is a zoom close to 0, showing data for A = —3/4 and 
— 1/2. The inset (b) shows the prefactor a vs A extracted 
from the numerical data (circles), for a larger set of values of 
A, which is compared with tv/2v = fi/ sin fi, fi = cos -1 A. 

Based on DMRG data obtained on critical open chains 
of size 200 < L < 1000, we find this proportionality is 
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still true even when A ^ where the decaying behavior 
of Ea is controlled by K. More precisely, plotting Sa as 
a function of — Ea for various values of the anisotropy A 
in Fig. 01 we find a linear relation Sa = — olEa with a 
prefactor perfectly described by a = /i/sin/i, as shown 
in the inset of Fig. with /i = cos -1 A. We note that 
the spin-wave velocity for the XXZ model is given by 
v = 7r(sin /i)/(2/Lt) so that we may write this relation as 



S A 



-(tto 2 /2v)E A , 



(13) 



where we have introduced the lattice spacing, a, to make 
the EE a dimensionless quantity (Ea has dimensions of 
energy per unit length.). We emphasize that Eq. I|13|) 
even holds for the Heiscnberg model with a = 1 where we 
find both Ea(L/2, L) and Sa(L/2, L) display the same 
logarithmic corrections, Eq. 10, as shown in Fig. 03 The 
sum Ea(L/2, L) + Sa(L/2, L) is found to rapidly decay 
as a power-law, with a power ~ 2.59 (see Fig. |5J. In- 
terestingly, for the SU(2) invariant model with J\ = 1 
and J2 — 0.241167, again linearity is observed, but with 
a prefactor a ~ 1.001689 not related to the spin veloc- 
ity, v, which we have determined to be v ~ 1.1699 [2fl |. 
Hence, Eq. I|13|) does not hold in this case. In Fig. [S] we 
also show the sum Ea(L/2, L) + aSA{L/2, L), for this 
model, which decays rapidly with a power ~ 2.56. 
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FIG. 5: Comparison between the alternating part of EE 
Sa and Ea [Eq. Q] from DMRG with m = 512 states, 
for Heisenberg models. Power-law decay of Sa(L/2, L) + 
oEa(L/2, L) drawn in a log-log plot, with a = 1 for the near- 
est neighbor chain (J2 = 0) and a = 1.00169 at the critical 
second neighbor coupling J2 = 0.241167. Lines are power-law 
fits: ~ IT 2 ' 56 for J 2 = 0.241167 and L -2 ' 59 for J 2 = 0. 



Conclusion — We emphasize that this alternating term 
in S(x, L) is universal and should not be regarded as a 
correction due to irrelevant operators. First of all, it is 
not a "correction" , since it is alternating. Secondly, it de- 
cays with the same power law as Ea{x,L) which is seen 



to be a property of the fixed point, not the irrelevant 
operators. (However, for the Heisenberg model, A = 1, 
the log factor in Ea(x,L) is due to the marginally irrel- 
evant operator.) The presence of a universal alternating 
term in S(x, L) is connected with the antiferromagnetic 
nature of the Hamiltonian (no t appearing, for example, 
in the quantum Ising chain |l0() and does not seem to fol- 
low from the general CFT treatment in Q ■ An analytic 
derivation of this phenomena remains an open problem. 
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